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A  USER’S  GUIDE  TO  THE  SEVP  SOLVER: 

AN  EFnCIENT  DIRECT  SOLVER  FOR  ELLIPTIC 
PARTIAL  DIFFERENTIAL  EQUATIONS 


1 .  INTRODUCTION 

This  report  describes  the  Stabilized  Error  Vector  Propagation  (SEVP) 
fortran  subroutines,  which  have  been  developed  at  the  Naval  Research 
Laboratory  (NRL)  for  the  solution  of  elliptic  partial  differential  equations 
using  finite  difference  methods.  The  SEVP  subroutines  have  been  used 
extensively  with  the  NRL  numerical  weather  prediction  models.  Since  the  code 
was  first  developed  (Madala,  1978),  various  enhancements  have  been  made.  The 
code  is  now  also  being  more  widely  used  at  NRL  and  at  other  research 
institutions  and  a  users  guide  has  therefore  become  necessary. 

The  SEVP  method  is  an  efficient  direct  method  which  is  used  to  solve  the 
finite  difference  approximation  to  an  elliptic  partial  differential  equation 
of  the  form 

a(x,y)  +  b(x,y)  +  c(x,y)  ^ 

8x^  0y  9x 

+  d(x,y)  ^  +  e(x,y)  ^  =  f(x,y) 

8y 

in  which  a(x,y)  b(x,y)  >0  on  a  rectangular  domain  in  some  arbitrary  x  and 
y  coordinate  system.  Grid  points  may  be  of  variable  spacing.  The  boundary 
conditions  for  ^(x,y)  can  be  Dirichlet,  Neumann,  periodic  in  one  direction, 
or  a  mixed  combination  of  these.  The  SEVP  solver  can  be  used  for  separable  and 
non-separable  elliptic  equations. 

Common  examples  of  these  elliptic  partial  differential  equations  are 
Poisson’s  equation  V2^ 

=  f  in  two  dimensions  and  the  special  case  when 
f  =  0,  Laplace’s  equation  in  two  dimensions.  For  the  boundary  conditions  in 
each  case,  the  value  of  ^  (Dirichlet  boundary  condition)  or  the  normal 
derivative  of  f  (Neumann  boundary  condition)  are  specified  on  each  boundary. 
Poisson’s  equation,  for  example,  arises  in  boundary  value  problems  in  which  a 
stream  function  is  used  to  solve  for  the  motion  of  an  incompressible  fluid 
when  the  vorticity  is  given.  A  Poisson’s  equation  is  also  encountered  whenever 
a  steady  state  potential  distribution  is  produced  in  response  to  sources  or 
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sinks  of  a  quantity,  where  the  rate  of  flow  of  the  quantity  is  proportional  to 
the  gradisnt  of  the  potential  distribution. 

As  shown  in  Madala  (1978),  the  SEVP  method,  which  is  based  on  the  error 
vector  propagation  (EVP)  method  (Hirota,  et  al.,  1970;  for  example),  is  much 
faster  than  successive  over-relaxation  (by  10  to  20  times)  and  uses  an  order 
of  magnitude  less  storage  space  than  the  method  of  Lindzen  and  Kuo  (1969). 
Unlike  the  EVP  method,  which  was  shown  by  McAvaney  and  Leslie  (1972)  to  be 
unstable  for  more  than  40  grid  points  in  the  marching  direction,  the  SEVP 
method  is  stable  for  any  number  ot  grid  points.  The  stabilization  is 
accomplished  by  dividing  the  integration  region  into  EVP  stable  overlapping 
blocks  and  imposing  artificial  boundaries  between  the  blocks. 

In  this  report,  the  form  of  the  difference  equation,  boundary  conditions 
and  the  subroutine  arguments  for  the  SEVP  solver  are  first  described.  We  then 
describe  the  derivation  of  the  difference  equation  and  boundary  conditions 
from  a  two-dimensional,  elliptic  partial  differential  equation,  defined  on  a 
rectangular  grid.  The  following  sections  describe  the  stabilized  error  vector 
propagation  method  and  each  of  the  SEVP  subroutines.  In  the  initial  set  up  in 
subroutine  SEVP,  the  difference  equations  are  modified  for  the  boundary 
conditions.  The  sections  on  the  preprocessor  subroutine  BSMl  and  the  solution 
subroutines  BSM2,  BSM3  may  be  skipped  by  the  reader  who  is  only  interested  in 
using  the  routines.  As  an  example,  Poisson’s  equation  is  solved  in  spherical 
coordinates  on  a  domain  on  the  surface  of  the  earth  bounded  by  constant 
latitude  and  longitude  circles.  Two  test  examples  demonstrating  the  use  of 
Dirichlet  and  Neumann  boundary  conditions  in  the  y  (marching)  direction  with 
periodic  boundary  conditions  in  x  are  described.  A  listing  of  the  test 
programs  is  provided  in  the  appendix. 
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2.  SUBROUTINE  SEVP 


Subroutine  SEVP  solves  a  system  of  linear  equations  of  the  form 

AX(i,j)*X(i-l,j)  +  Ay(i,j)*X(i.j-l)  +  BB(i, j )*X(i. j ) 

+  CX(i,  j)*X(i+l.j)  +  CY(i,  j)*X(i,  j+1)  =*  F(i,j)  (2.1) 

for  i  =  2,...,M-1  and  j  =  2,...,N-1. 

with  boundary  conditions  for  j=l,  j=N,  i=l  and  i  =  M  of  the  form 
X(i,l)  =  (l-All)*X(i,l)  +  All*(X(i,2)  -  Fll(i)) 

X(i,N)  =  (l-AlN)*X(i,N)  +  AlN*(X(i,N-l)  +  FlN(i))  for  i  =  1 . M 

and 

X(l,j)  =  (1-A21)*X(1, j)  +  A21*(X{2.j)  -  F21(j))  (2.2) 

X(M,j)  -  (1-A2M)*X(M. j)  +  A2M*(X(M-l,j)  +  F2M(j))  for  j  -  1 . N 

or  X(l,j)  *  X(M-l,j) 

X(M,j)  =  X(2.j)  for  j  =  1.....N. 

These  difference  equations  result  from  the  discretization  of  elliptic 
partial  differential  equations  of  the  form 

a(x,y)  +  b(x,y)  +  c(x,y)  ^ 

9x  9y  9x 

+  d(x,y)  ^  +  e(x,y)  ^  *  f(x,y)  (2.3) 

9y 

defined  in  a  rectangular  domain  x^  <  x  <  x^,  Fa  X  Xb  which 

a(x,y)  b(x,y)  >  0,  and  the  x  and  y  represent  some  coordinate  system.  Grid 
points  can  be  of  variable  spacing.  Boundary  conditions  for  ^(x,y)  may  be 
Dirichlet  or  Neumann  at  the  x  and  y  boundaries  x  =  Xg,  x  =  X5, 
y  =  yg,  y  =  y^;  or  periodic  in  x;  or  a  mixed  combination  of  these. 
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2.1  Description  of  Arguments 

The  SEVP  subroutines  are  called  by  the  statement 

CALL  SEVP ( AX , AY , CX , CY, BB , RINV , RINVl , RCOR, NBSIZ2 , IS , IE , INX , SUMF , 
DOMMYl ,  RTILDA ,  X ,  F ,  ERR ,  FI  1 ,  FIN ,  F2 1 ,  F2M ,  M,  M .  Ml ,  N1 ,  M2 ,  N2 ,  NBLK , 

NBLKl , B1 , B2 , All , AIN , A21 , A2M. IC Y , TOL , IFLG ) 

where  the  arguments  are  as  follows: 

QR.  iiffqi 

M 

The  number  of  points  in  the  i  direction,  including  the  boundary. 
Ml  =  M-1,  M2  =  M-2. 


The  number  of  points  in  the  j  (marching)  direction,  including  the  boundary. 
N  must  be  larger  than  4. 

N1  =  N-1,  N2  =  N-2. 

NBLK 

Represents  the  number  of  blocks  in  j  (marching)  direction.  For  N  larger 
than  14,  NBLK  =  (N+3)/9  gives  blocks  of  approximately  8  to  13  rows  each. 
For  N  less  than  11  but  larger  than  4,  use  NBLK  =  1;  while  for  N  between  11 
and  14,  one  or  two  blocks  may  be  used. 

Note:  For  the  case  of  Neumann  boundary  conditions  at  boundary  j  ®  N  (i.e. 
AIN  =  1)  and  with  NBLK  larger  than  one,  last  block  is  fixed  with  7  rows  by 
SEVP  solver,  to  allow  for  a  more  efficient  convergence  of  the  solution.  In 
this  case,  use  formula  to  calculate  the  number  of  blocks  required  for  the 
first  N-5  rows  and  then  add  one  to  get  total  number  of  blocks. 

NBLKl 

=  NBLK-1,  for  NBLK  greater  than  one. 

=  1  for  NBLK  equal  to  one. 

AX,  AY,  CX,  CY,  BB 

Two-dimensional  arrays  of  dimension  M  x  N  containing  the  coefficients  of 
the  difference  equation. 


X 

Two-dimensional  array  that  on  input  contains  an  initial  guess  for  solution 
on  interior  points  i  =  2,... M-1,  j  =  2,..., N-1;  and  for  Dirichlet  boundary 
conditions,  contains  the  boundary  values  at  the  appropriate  boundaries. 
Initial  guess  might  be  an  average  of  the  boundary  values  (for  Dirichlet 
BC’s)  or  a  zero  value  (for  Neumann  BC’s),  for  example. 


Two-dimensional  array  of  dimension  M  x  N  that  specifies  the  forcing  on  the 
right  hand  side  of  the  difference  equation. 
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Fll,  FIN 

One- dimensional  arrays  of  length  M  that  contain  coefficients  for  the 

Neumann  boundary  conditions  at  the  j  boundaries. 

F21,  F2M 

One-dimensional  arrays  of  length  N  that  contain  coefficients  for  the 

Neumann  boundary  conditions  at  the  i  boundaries. 

All,  AIN,  A21,  A2M 

Real  variables  which  take  the  value  0.0  for  Dirichlet  boundary  conditions, 
and  1.0  for  Neumann  boundary  conditions  at  their  respective  boundaries.  All 
corresponds  to  j  =1,  AIN  to  j  =  N,  A21  to  i  =  1,  and  A2M  to  1  «  M.  For 

periodic  boundary  conditions  in  i  (with  ICY  =  1),  A21  and  A2M  are  zero. 


ICY 

=  1  for  periodic  boundary  conditions  in  i. 
=  0,  otherwise. 


TOL 

Error  tolerance  for  SEVP  solver. 

IFLG 

=  0  for  first  call  of  SEVP  solver.  Preprocessor,  subroutine  BSMl  called. 

=  1  for  subsequent  calls  if  coefficients  AX,  AY,  CX,  CY,  BB  of  difference 
equation  and  the  boundary  condition  type  All,  AIN,  A21,  A2M,  ICY  do  not 
change.  Preprocessor  is  then  skipped. 

ON  OUTPUT 

X 

Two-dimensional  array  of  dimension  M  x  N  containing  the  solution. 

ERR 

Two-dimensional  array  of  dimension  M  x  N  containing  the  residual  error  for 
the  difference  equation. 

OUTPUT  ARRAYS  COMPUTED  BY  PREPROCESSOR 

These  arrays  must  not  be  destroyed  if  SEVP  will  be  called  again  with 
IFLG  =  1. 

RINV 

Three-dimensional  array  of  dimension  M2  x  M2  x  NBLK,  containing  two- 
dimensional  M2  X  M2  influence  or  residual  matrices  for  the  blocks.  The 
first  NBLK-1  influence  matrices  relate  the  solution  error  on  the  last 
interior  row  of  each  block  to  that  on  the  second  row  of  that  block.  The 
last  two-dimensional  matrix  is  a  residual  matrix  relating  the  residual 
error  computed  on  the  last  interior  row  to  the  solution  error  on  the  second 
row  of  the  last  block. 

RINVl 

Three-dimensional  array  of  dimension  M2  x  M2  x  NBLKl,  containing  NBLK-1 
two-dimensional  M2  x  M2  influence  matrices,  which  relate  the  solution  error 
on  the  last  two  rows  of  each  block,  except  the  last. 
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WORKING  INTEGER  ARRAYS 
NBSIZ2 

One-dimensional  integer  array  of  size  NBLK,  containing  the  number  of 
interior  rows  in  the  j  (marching)  direction  between  the  end  boundaries  of 
the  blocks.  The  first  block  contains  a  total  of  NBSIZ2(1)  +  2  rows,  while 
the  remaining  blocks  contain  NBSIZ2()  +  3  rows  each. 

IS,  IE 

One-dimensional  integer  arrays  of  size  NBLK,  defining  the  beginning  and  end 
row  of  each  block,  respectively. 


INX 

One -dimensional  integer  array  of  size  M  containing  i  indices. 

WORKING  REAL  ARRAYS 
SUMF 

One-dimensional  array  of  size  NBLK  containing  the  sum  of  the  absolute 
values  of  the  forcing  on  the  last  interior  row  of  each  block. 

RCOR 

Two-dimensional  array  of  dimension  M  x  3  containing  set  of  three  row 
vectors  used  for  marching  in  j  direction. 

RTILDA 

One-dimensional  error  vector  of  size  M-2. 

Bl,  B2 

One-dimensional  work  arrays  of  size  M-2. 

DUMMYl 

Two-dimensional  work  arrays  of  size  M-2  x  M-2. 

INTERNAL  PARAMETERS  IN  SUBROUTINE  SEVP 

lOUT 

To  print  residual  errors  after  each  iteration  set  iout=l,  otherwise  iout=0. 

EPS 

The  machine  precision  of  computer,  used  as  a  lower  limit  on  the  tolerance 
TOL.  For  the  Cray  with  64  bit  words,  single  precision  floating  point  words 
(with  48  bit  mantissa)  are  precise  to  14  decimal  places  or  for  a  machine 
precision  of  about  0.5  x  10"^^. 

MAXSIZ 

Maximum  number  of  points  allowed  in  each  block.  For  the  precision  of  the 
Cray  computer  setting  MAXSIZ  =  15  ensures  the  solver  converges  rapidly  to  a 
solution. 
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IX,  JX.  JF 

For  4-8lded  Neuaiann,  or  Z-slded  Neumann  with  cyclic  boundary  conditions,  a 
boundary  value  must  be  prescribed  at  point  (IX,  JX)  for  a  solution  to  be 
obtained.  JF  is  the  nearest  interior  row  to  the  prescribed  boundary  point 
(IX,  JX).  For  boundary  point  IX-2,  JX-1,  the  nearest  interior  row  is  JF-2. 

PROGRAM  LANGUAGE  Fortran-77 

PRECISION  Single  precision  on  64  bit  computer.  For  32  bit  computer, 

subroutines  should  be  modified  for  double  precision. 

ORIGINATOR  Rangarao  V.  Madala 

Naval  Research  Laboratory,  Washington  D.C.  20375 


7 


3 .  DIFFERENCE  APPROXIMATION 


In  this  section  we  show  how  the  difference  equation  is  derived  from  a  two- 
dimensional  elliptic  partial  differential  equation.  A  two-dimensional  elliptic 
equation  of  the  form 


a(x.y)  ^ 
dx^ 

+  b(x,y)  ^ 

+  d(x,y)  ^ 

8y 

8^ 

+  c(x,y) 

8x 

+  e(x,y)  ^ 

f(x,y) 

(3.1) 

is  assumed 

to  be  defined 

in  a  .ectangular 

domain 

Xg  <  X  <  xb,  and 

Yjj  <  y  <  y^j,  where  a(x,y)  c(x,y)  >  0  ,  with  boundary  conditions  specified  at 
boundaries  x=Xa.  x=Xb*  y~7a'  y=yb*  Boundary  conditions  may  be  Dirichlet, 
Neumann  in  x  and  y  or  periodic  in  x. 


As  an  illustration,  we  define  a  grid  mesh  of  M  by  N  points  (x^,  y j )  with  a 
grid  spacing  of  constant  Ax  and  Ay.  However,  variable  grid  spacing  can  be 
used,  resulting  in  modification  to  the  coefficients  of  the  difference 
equation.  A  second  order  finite  difference  approximation  of  Eq.  (3.1),  using  a 
five  point  stencil,  is  then 


a .  . 
i.J 


4i 


i+li 


+  b.  ^^.Hl  '  ^  ^  i,j _ ^ 

Ay^ 


+  c,  .  — 

i.J 

2  Ax 

+  a .  . 

i.J 

+  e 

i.j 

= 

where 

^i.j  ■ 

^(x.  , 

Yj). 

It 

f  (x^. 

and 

a .  .  = 

i.J 

a(x^. 

y^),  etc 

By  defining 

X(i,  j) 

and 

2  Ay 


f .  . 
i.J 


(3.2) 
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AX(i,j) 


AY(i.j)  = 

CX(i,j)  = 

CY(i,j)  = 

= 

F(i.j)  = 


Ax^ 
b 


LI 


Ay 

a 


LI 


Ax 

b 


ill 


Ay^ 

2  a.  . 
_ LlI 

Ax^ 


lia 

2  Ax 

■ 

2  Ay 
c .  . 

2  Ax 

2  Ay 

2  b.  . 
_ LLL 

Ay^ 


+  e . 


i.J 


'i.  j 


Our  difference  equation  (3.2)  can  then  be  written  in  the  form 
AX(i.j)*X(i.l.j)  +  AY(i,j)*X(i.j-l)  +  BB(i.j)*X(i.j) 

+  CX(i,j)*X(i+l.j)  +  CY(i,j)*X(i,j+l)  -  F(i,j) 
for  i  =  2 . M-1  and  j  =  2,...,N-1. 


(3.3) 


(3.4) 


We  will  now  demonstrate  the  use  of  Neumann,  Dirichlet  and  periodic 
boundary  conditions  in  x. 


3 . 1  Dirichlet  Boundary  Conditions 

For  Dirichlet  boundary  conditions  in  x,  boundary  values  are  given  at  and 

X=Xb 

/'(xa-y)  “  p(y) 

and  ^(Xb*y)  =  'l(y)  for  Ya  ^  Y  ^  /b 

For  our  grid  mesh,  defining  M  points  in  x  with 
Xi  =  (i-1)  Ax  +  Xg  for  i  =  1,...,M 
where  Ax  =  (x^-Xg) / (M-1) ,  we  have 
-  Pj 

^(M,j)  »  qj  for  j  =  1 . N 

By  defining  A21  =  A2M  =  0.0  and 

X(l,j)  =  Pj  and  X(M,j)  =  qj  for  j=l,...,N 
our  boundary  conditions  are  in  the  appropriate  form  of  Eq.  (2.2). 
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3 . 2  Neumann  Boundary  Conditions 

For  Neumann  boundary  conditions  in  x,  derivatives  are  specified  at  the 
boundaries  x=x^  and  x=xij 

—  (X  .y)  =  u(y)  and  —  =  v(y)  for  y  <  y  <  y. 

3x  ®  8x  ^  a  b 


For  our  grid  mesh,  we  include  two  columns  of  fictitious  points  at  x^-Ax/Z  and 
Xb+Ax/2,  and  define  our  x  grid  points  with  a  shift  of  Ax/2,  so  that 

xi  =  (i-1)  Ax  +  Xg  -  Ax/2  for  i  =  1, _ ,M 

with  Ax  =  (x^-Xa) / (M-2) .  Then  using  centered  differences  we  have 


u .  and 
1 


Ax 


V.  for  j  =  2 . N-1 

J 


By  defining  A21  =  A2M  =  1.0, 

F21(j)  =  Ax  Uj ,  and  F2M(j)  =  Ax  Vj  for  j  =  2,... N-1 
we  have  our  boundary  conditions  in  the  appropriate  form  of  Eq.  (2.2). 


3 . 3  Periodic  Boundary  Conditions 

Suppose  that  ^  is  a  periodic  function  in  x  with  period  L  =  x^-Xg  so  that 

^(x+xb,y)  =  ^(x+Xa,y) 

With  a  column  of  points  included  at  x=X5+Ax,  we  define  our  grid  mesh  over  a  x 
domain  of  Xg  ^  x  ^  x^+Ax  so  that 
x^  =  (i-1)  Ax  +  Xg  for  i  =  1,...,M 
with  Ax  =  (x^-Xg) / (M-2) .  Then  we  have 


M-1 


j 


and 


"M.j 


2.j 


for  j  =  1, . . . ,N. 


Setting  ICY  =  1  with  A21  =  A2M  =  0.0  will  give  the  boundary  conditions 
X(l,j)  =  X(M-l,j) 

X(M,j)  =  X(2,j),  for  j  =  1,...,N 
as  required  in  Eq.  (2.2). 
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4.  SOLUTION  METHOD 


The  method  is  based  on  the  error  propagation  method  (EVP)  or  sweep  out 
method  (Hirota  et  al.,  1970).  In  the  EVP  method,  a  forward  sweep  provides  a 
particular  solution,  which  is  followed  by  a  further  forward  sweep  of  the 
solution  error  and  the  correction  of  the  solution.  We  demonstrate  the  method 
for  the  case  of  Dirichlet  boundary  conditions. 

At  the  interior  points,  we  start  with  an  initial  guess  solution,  which  may 
be  a  constant  value  equal  to  the  average  of  the  boundary  values,  for  example. 
Rearranging  the  difference  equation  (2.1)  we  can  write 

X(i,j+1)  =  [F(i.j)  -  AX(i,j)*X(i-l,j)  -  AY(i,j)*X(i, j-1)  -  BB(i.j)*X(i.j) 

-  CX(i,j)*X(i+l,j)]  /  CY(i.j) 

for  i  =  2 . M-1  and  j  =  2 . N-1.  (4.1) 

A  particular  solution  Xp(i,j)  can  be  found  by  starting  with  the  initial  guess 
on  the  second  row  j  =  2  and  marching  Eq.  (4.1)  forward  row  by  row  to  the  end 
boundary  j  =  N.  The  particular  solution  Xp(i,j)  then  satisfies  the  boundary 
conditions  on  three  sides  j=l,  i=l,  i=M  but  differs  from  the  solution  at 
the  end  boundary  j  =  N  by  an  error  row  vector  defined  by 
CN(i)  =  X(i,N)  -  Xp(i,N) 

=  X(i,N)  -  [F(i,N-l)  -  AX(i,N-l)*Xp(i-l,N-l)  -  AY(i,N-l)*Xp(i,N-2) 

-  BB(i,N-l)*Xp(i,N-l)  -  CX(i,N-l)*Xp(i+l,N-l) ]  /  CY(i,N-l) 
for  i  =  2 . M-1.  (4.2) 

A  solution  error  Xh(i.j)  =  X(i,j)  -  Xp(i,j)  can  be  defined  which  satisfies 
the  homogeneous  difference  equation  (Eq.  (2.1)  with  F  =  0)  with  the 
homogeneous  boundary  conditions  on  the  three  sides 
Xh(i,l)  =  0,  for  i  =  1,...,M 

Xhd.j)  =  Xh(M,j)  =0  for  j  =  1 . N  (4.3) 

and  a  non-zero  boundary  value  given  by  the  solution  error 

Xh(i,N)  =  ^  =  2,...,M-1  (4.4) 

at  the  end  boundary  j  =  N.  We  can  compute  an  influance  matrix  which  relates 
the  solution  error  on  the  last  interior  row  to  that  on  the  second  row.  The 
solution  error  can  then  be  obtained  for  the  whole  domain  by  marching  forward 
from  the  second  row  and  the  particular  solution  corrected  to  give  thj^ 
solution. 

To  obtain  the  influence  matrix  we  start  with  a  unit  row  vector  ej^  defined 
on  the  second  row.  The  unit  row  vector  e^  has  the  element  unity  and  the 
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remaining  elements  zero.  Marching  the  unit  row  vector  forward,  using  the 

homogeneous  difference  equation,  we  obtain  influence  row  vectors  tj^  on  the 
last  row  j  =  N,  defining  the  row  of  our  influence  matrix  T.  Repeating  the 
marching  sequence  for  all  the  M-1  unit  row  vectors  defines  a  set  of  M-1 
influence  row  vectors,  completing  the  influence  matrix  T.  Using  a  linear 

combination  of  influence  vectors  and  the  principle  of  linear  superposition,  we 

can  relate  our  error  row  vector  (for  the  last  row  j  =  N)  to  the  error 

vector  “  Xh(i.2)  )  on  the  first  interior  row  j  =  2.  That  is,  in  matrix 

form  using  the  influence  matrix,  we  have 

eN  =  €2  T  (4.5) 

Therefore  given  the  computed  particular  solution  Xp(i,j)  and  the  error 
vector  Cn  on  the  last  row,  the  error  vector  can  be  computed  on  the  first 
interior  row, 

€2  =  Cn  T-1.  (4.6) 

Then  marching  the  solution  error  forward  provides  the  solution  error 
throughout  the  whole  domain,  and,  when  combined  with  the  particular  solution, 
completes  the  solution. 

The  niarching  of  the  solution  error  can  alternatively  be  conceptualized  as 
follows.  The  error  vector  ejj  at  the  end  boundary  j  »=  N  can  be  related  to  a 

residual  vector  p  defined  on  the  last  interior  row  j=N-l.  Rearranging  Eq. 

(4.2),  we  have 

AX(i,N-l)*Xp(i-l,N-l)  +  AY(i,N-l)*Xp(i,N-2)  +  BB(i,N-l)*Xp(i,N-l) 

+  CX(i,N-l)*Xp(i+l,N-l)  +  CY(i,N-l)*X(i,N) 

=  F(i,N-l)  +  CY(i,N-l)*eN(i)  (4.7) 

where  we  define  the  residual  vector  p  on  the  last  interior  row  j  =  N-1  by 

/J(i)  =  CY(i,N-l)*eN(i) ,  for  i  =  2,...,M-l.  (4.8) 

The  solution  error  can  now  be  thought  of  satisfying  the  homogeneous  boundary 
conditions  on  all  sides  and  the  homogeneous  difference  equation  at  all  rows 
except  at  the  last  interior  row.  At  the  last  row  j=N-l,  the  solution  satisfies 
the  difference  equation  with  a  forcing  function  given  by  the  residual  p.  By 
similarly  marching  the  unit  vectors  forward,  a  set  of  residuals  can  be 
computed  on  the  last  interior  row,  defining  a  residual  matrix  R  (Hirota,  et 
al,  1970).  The  residual  matrix  then  relates  the  solution  error  on  the  second 
row  to  the  residual  p  defined  on  the  last  interior  row. 

^  =  62  R  («.9) 
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The  two  equivalent  views  of  the  solution  error  arc  shown  schematically  below 
in  Fig.  1. 

Due  to  round-off  errors  in  marching  the  particular  solution  forward  and 
due  to  inaccuracies  in  the  computed  inverse  T“^,  the  error  vector  62  computed 
on  the  second  row  using  Eq.  (4.6)  is  an  approximation  to  the  true  error 
vector.  Marching  the  error  vector  forward  then  and  correcting  the  particular 
solution  provides  an  approximate  solution.  Provided  that  the  error  vector  is 
not  too  inaccurate,  an  improved  solution  can  be  found  by  recomputing  a  new 
error  vector  on  the  last  row  and  repeating  the  above  procedure.  A  number  of 
iterations  of  the  procedure  are  usually  required  to  converge  to  an  accurate 
solution. 


SOLUTION  ERROR 


i=l  6=0  i=M  e  *  0  i-M 


- ? - 

1 

1 

1 

1 

1 

1 

1 

1  « 
tv> 

1 

1 

1 

1 

<®k>  Cr 

H 

0 

0 

T  =  {tfc} 

R  -{rfe}.  p  -  €2  R 

e  =  €n  =  C2  t 

e  -  0 

Fig.  1.  The  two  equivalent  views  of  the  boundary  conditions  and  forcing  for 
the  solution  error,  (a)  Solution  error  satisfies  homogeneous  difference 
equation  with  a  non  zero  boundary  condition  at  j=N  given  by  the  error  vector 
Cn-  Influence  matrix  T  found  by  marching  set  of  unit  vectors  {ej^}  to  last  row 
j  *  N.  (b)  Solution  error  satisfies  homogeneous  boundary  conditions  with  a 
non-zero  line  forcing  function  given  by  the  residual  p  at  j=N-l.  Residual 
matrix  R  found  by  marching  unit  vectors  {e^}  to  last  interior  row  j  =  N-1  and 
computing  residuals  on  that  row. 
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As  pointed  out  by  McAvaney  and  Leslie  (1972),  the  residual  matrix  (or 
equivalently  the  Influence  matrix)  is  ill-conditioned,  with  the  degree  of  ill- 
conditioning  increasing  rapidly  with  increasing  number  of  points  N  in  the 
marching  direction.  For  large  N,  the  inverse  of  the  influence  matrix  may  not 
be  computed  very  accurately  and  for  moderate  round-off  error  any  accuracy  in 
the  computed  error  vector  may  be  lost.  The  iterative  procedure  then  will  not 
converge  to  a  solution.  The  accuracy  of  the  inverse  therefore  limits  the 
number  of  points  that  can  be  effectively  used  in  the  marching  and  correction 
process.  On  the  Cray  computer  with  64  bit  words,  this  limit  is  about  20 
points.  The  degree  of  ill-conditioning  of  the  influence  matrix,  which  leads  to 
inaccuracies  in  the  inverse,  can  be  measured  by  a  condition  number  defined  by 
cond(T)  =  I iTl I  I IT-^I I  (4.10) 

where  llTlI  and  llT“^ll  are  matrix  norms  of  the  influence  matrix  and  it’s 
inverse,  respectively  (see  Dahlquist  and  Bjorck,  1974;  for  example).  Then 
given  an  error  llffe({ll  in  the  error  vector  on  the  last  row,  an  upper  bound  on 
the  relative  error  of  the  computed  error  vector  on  the  second  row  can  shown  to 
be 


\\6eJ\ 

-  i  cond(T)  -  (4.11) 

lle^ll  116^11 

where  the  matrix  and  vector  norms  are  defined  in  a  consistent  manner  (see 
Dahlquist  and  Bjorck).  For  example,  taking  the  case  of  a  domain  with  14  points 
in  the  marching  j  direction  and  51  points  in  i  direction  the  condition  number 
is  computed  to  be  of  order  10®.  Then  given  an  extreme  case  of  a  high  relative 
round-off  error  of  10“^2  in  the  computed  error  vector  Cn  on  the  last  row,  when 
using  single  precision  floating  point  numbers  on  the  64  bit  Cray  computer,  we 
expect  the  computed  error  vector  on  the  second  row  to  have  a  precision  of  at 
least  four  significant  figures  (for  a  relative  error  of  10“^).  The  iterative 
procedure  will  then  converge  to  a  solution. 

To  stabilize  the  EVP  method  for  a  large  number  of  grid  points,  Madala 
(1978)  divided  the  domain  into  a  series  of  blocks  in  the  marching  direction. 
Adjacent  blocks  overlap  in  the  marching  direction  so  that  the  first  and  last 
two  rows  of  adjacent  blocks  are  common.  To  provide  for  a  rapid  convergence 
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i-1  e  -  0  i-M 


Fig.  2.  The  relationship  between  the  error  vectors  6  on  various  rows  at  the 
block  boundaries  and  the  residual  error  ^3  on  the  next  to  the  last  row  jaN-1, 
for  the  case  of  three  blocks.  Rn,  R12*  *^21*  ^22  influence  matrices  and 
R3]^  Is  a  residual  matrix,  defined  on  the  rows  shown. 


of  the  iterative  procedure,  the  number  of  points  in  each  block  is  limited  to 
less  than  15  points. 

In  the  preprocessor  step,  NBLK  x  2  -  1  matrices  are  computed  which  relate 
the  solution  error  on  various  boundary  rows  of  each  block.  In  Pig.  2,  showing 
three  blocks  for  example,  the  influence  matrices  Ri2t  Rn  for  the  first  block 
relate  the  error  vector  €1  on  the  second  row  to  the  error  vector  €3  on  the 
last  row  and  to  the  error  vector  €2  on  the  next  to  the  last  row,  respectively. 
Similar  influence  matrices  R22*  *^21  computed  for  the  second  block.  For  the 
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last  block  a  residual  matrix  R31  relates  the  residual  solution  error  ^3  on  the 
last  interior  row  j=N-l  to  the  error  vector  £5  on  the  second  row  of  the  last 
block. 

In  the  forward  sweep  a  particular  solution  is  found  by  the  marching  and 
correction  process  above,  satisfying  the  boundary  conditions  and  the  initial 
guess  field  on  the  interior  block  boundaries.  The  particular  solution 
satisfies  the  forcing  function  at  all  points  except  on  the  last  interior  row, 
where  the  difference  equation  satisfies  the  forcing  plus  a  residual  ^3.  In  the 
backward  sweep  a  solution  error  can  be  computed  on  the  block  boundaries  and 
the  first  interior  row  from  the  residual  on  the  last  interior  row,  using 
the  residual  and  influence  matrices  defined  for  the  blocks.  By  marching  the 
solution  error  forward  in  each  block  the  particular  solution  is  then  corrected 
to  give  the  solution.  The  details  of  the  stabilized  error  vector  propagation 
method  are  described  in  the  next  chapter. 
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5.  DESCRIPTION  OF  SEVP  SUBROUTINES 


5.1  Initial  Set-up.  SEVP 

In  the  main  subroutine  SEVP,  the  block  boundaries  are  first  computed  for 
the  prescribed  number  of  blocks  and  the  difference  equations  are  modified  to 
accommodate  boundary  conditions  which  are  not  Dirichlet. 

Fig.  2  shows  the  relationship  of  the  block  boundaries  for  the  case  of 
three  overlapping  blocks.  The  beginning  and  end  rows  of  the  blocks  are  defined 
by  the  one -dimensional  arrays  IS  and  IE,  respectively.  The  integer  array 
NBSIZ2,  defining  the  number  of  interior  rows  between  the  end  boundaries  of  the 
blocks,  is  computed  for  the  prescribed  number  of  blocks  in  subroutine  BLKSIZ. 
The  first  block  contains  a  total  of  NBSIZ2(1)  +  2  rows,  while  the  remaining 
blocks  contain  NBSIZ2()  +  3  rows  each.  For  the  special  case  of  Neumann 
boundary  conditions  at  boundary  j  =  N  (i.e.  AIN  ■  1),  the  last  block  (third 
block  in  Fig.  2)  is  set  with  NBSIZ2()  =  4  by  BLKSIZ  (for  a  total  of  7  rows  for 
the  block),  to  provide  for  a  more  efficient  convergence  of  the  solution. 

For  the  case  of  cyclic  boundary  conditions  we  refine  an  i  index  INX(i)  by 

INX(i)  =  i  for  i=2,...,M-l 

INX(l)  =  M-1  (5.1) 

INX(M)  =  2 

The  difference  equation  (2.1)  is  then  modified 

AX(i,j)*X(INX(i-l),j)  +  AY(i,j)*X(i,j-l)  +  BB(i,j)*X(i,j) 

+  CX(i,j)*X(INX(i+l),j)  +  CY(i,j)*X(i,j+l)  »  F(i,j)  (5.2) 

for  i  =  2,..., M-1  and  j  =  2,...,N-1. 
to  provide  for  the  cyclic  boundary  conditions  at  the  boundaries  i  =  1  and 
i  =  M, 

X(l,j)  =  X(M-l,j) 

X(M,j)  =  X(2,j),  for  j  =  2,...,N-..  (5.3) 

In  the  case  of  Neumann  boundary  conditions  at  the  x  boundaries  i  =  1  and 
i  =  M,  our  boundary  conditions  are 

X(l,j)  =  X(2,j)  -  F21(j)  and 

X(M,j)  =  X(M-l,j)  +  F2M(j)  for  j  =  2 . N-1.  (5.4) 


Now  the  dif'^  rence  equation  at  1  =  2  Is 

AX(2,j)*X(l.j)  +  AY(2,j)*X(2,j-l)  +  BB(2, j ) *X(2, j )  +  CX(2, j )*X(3, j ) 

+  CY(2,j)*X(2,j+l)  =  F(2,j)  (5.5) 

for  j  =  2, . . . ,N-1. 

Substituting  our  boundary  condition  Eq  (5.4)  for  X(l,j)  in  Eq.  (5.5)  we  have 
AY(2,j)*X(2.j-l)  +  [BB(2,j)  +  AX(2, j ) ]*X(2, j )  +  CX(2, j )*X(3. j ) 

+  CY(2,j)*X(2,j+l)  =  [F(2,j)  +  AX(2,j)*F21(j)]  (5.6) 

for  j  =  2, . . . ,N-1. 

By  redefining  our  variables  such  that 
X(l.j)  0.0 

BB(2.j)  -»  BB(2,j)  +  AX(2,j)  (5.7) 

F(2,j)  F(2,j)  +  AX(2,j)*F21(j)  for  j  =  2 . N-1 

we  can  leave  the  form  of  the  difference  equation  the  same  as  for  the  case  of 


Dirichlet  boundary  conditions.  Similarly,  by  redefining 
X(M,j)  — »  0.0 

BB(M-l.j)  BB(M-l,j)  +  CX(M-l,j)  (5.8) 

F(M-l,j)  — ►  F(M-l,j)  -  CX(M-1, j)*F2M(j)  for  j  =  2 . N-1 


the  form  of  the  difference  equation  is  maintained  for  i=M-l.  We  have 
essentially  folded  our  boundary  conditions  into  the  difference  equation, 
ostensibly  replacing  Neumann  boundary  conditions  with  Dirichlet  boundary 
conditions  with  zero  boundary  values  X(l,j)  =  X(M,j;  =  0.0  at  the  boundaries 
i  =  1  and  i  =  M,  respectively.  Similar  modifications  are  made  for  Neumann 
boundary  conditions  at  the  y  boundaries,  j  =  1  and  j  =  N. 

In  the  case  of  4- sided  Neumann  boundary  conditions  or  2- sided  Neumann  with 
cyclic  boundary  conditions  in  x,  the  solution  is  indeterminate.  For  if  ^  is  a 
solution,  then  so  is  ^  +  K,  where  A  is  an  arbitrary  constant.  For  the  solver 
to  obtain  a  solution,  a  boundary  value  must  be  prescribed  at  some  chosen 
boundary  point  (IX,  JX) .  For  the  most  efficient  convergence  of  the  solution, 
it  is  recommended  that  a  zero  value  be  used  in  this  case  for  the  initial  guess 
and  for  the  constant  value  at  the  prescribed  boundary  point.  For  the  nearest 
interior  row  j  =  JF  to  the  prescribed  boundary  point,  the  boundary  point 
X(IX,JX)  and  the  coefficients  BB(IX,JF)  and  the  forcing  F(IX,JF)  in  the 
difference  equation  are  not  modified  by  the  SEVP  subroutine.  The  boundary 
point  initially  chosen  in  SEVP  is  IX  =2,  JX  =  1,  with  the  nearest  interior 
row  being  JF  =  2. 
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For  convenience,  we  further  modify  the  difference  equation  for  j*N-l  to 
include  the  Dirichlet  boundary  condition  at  j*N  in  the  forcing  function.  In 
this  case  the  value  X(i,N}  of  the  solution  variable  is  prescribed  on  the 
boundary  j=N.  The  difference  equation  for  j=N-l  is 

AX(i,N-l)*X(i-l,N-l)  +  AY(i,N-l)*X(i,N-2)  +  BB(i,N-l)*X(i,N-l) 

+  CX(i,N-l)*X(i+l,N-l)  +  CY(i,N-l)*X(i,N)  =  F(i,K-l)  (5.9) 

for  i  =  2, . . . ,M-1, 
which  we  rearrange  as 

AX(i,N-l)*X(i-l,N-l)  +  AY(i,N-l)*X(i,N-2)  +  BB{i,N-l)*X(i,N-l) 

+  CX(i,M-l)*X(i+l.N-l)  =  F(i.N-l)  -  CY(i,N-l) *X(i,N)  (5.10) 

for  i  =  2, . . . ,M-1 . 

By  redefining  our  forcing  function  at  j=N-l  such  that 

F(i,N-l)  — ►  F(i.N-l)  -  CY(i,N-l)*X(i,N)  for  i  =  2 . M-1  (5.11) 

we  can  assume  for  convenience  that  X(i,N)  =  0.0  at  the  end  boundary  j  =  N. 

5.2  The  Preprocessor.  BSMl 

In  the  preprocessor,  subroutine  BSMl,  NBLK  x  2  -  1  matrices  are  computed 

which  relate  the  solution  error  on  various  boundary  rows  of  each  block  and  the 
residual  error  on  the  last  interior  row  j  «  N-1.  The  solution  error  satisfies 
homogeneous  boundary  conditions  on  all  sides,  and  the  homogeneous  difference 
equation  at  all  interior  points  except  the  last  interior  row  j  »  N-1.  On  the 
last  Interior  row  the  solution  error  is  forced  by  a  residual. 

To  obtain  the  influence  matrices  for  the  first  block,  the  unit  vector 
(with  the  k^^  element  one)  defined  on  the  second  row  j  =  2,  is  marched  forward 
with  homogeneous  boundary  conditions  to  compute  row  vectors  on  the  last  two 
rows  of  the  block,  j  =  IE(1)-1,  j  =  IE(1).  These  define  the  k^^^  rows  for  two 
influence  matrices  (shown  as  and  R12  in  Fig.  2)  for  the  last  two  rows. 
Repeating  the  marching  procedure  for  the  other  unit  vectors  completes  the 
influence  matrices  R^  and  Ri2t  which  then  relate  the  solution  error  vector  Ci 
on  the  second  row  to  the  error  vectors  €2,  63  on  the  last  two  rows, 
respectively,  so  that 

€2  =  €1  Rii  (5.12) 

and 
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*3  =  Cl  Ri2  (5.13) 

By  matrix  multiplication  of  th«i  inverse  of  times  R2^2  ®  further  influence 

matrix  can  be  obtained,  which  relates  the  solution  error  vector  on  the  next 
to  the  last  and  the  last  rows  of  the  first  block.  For  by  the  elimination  of 
in  Eqs.  (5.12)  and  (5.13),  we  have 

€3  =  C2  Rii”^R12  “  Si-  (5.14) 

The  inverse  of  the  influence  matrix  is  computed  using  =  Ri2"^  ®^11* 

For  the  second  block  (and  any  other  subsequent  interior  blocks),  we  again 
start  with  a  unit  vector  defined  on  the  second  row  of  the  block,  but  the 
vector  on  the  first  row  is  now  defined  by  the  k^^  column  of  the  inverse  of  the 
influence  matrix  defined  for  the  previous  block.  For  if  ej^  is  the  unit 
vector,  then  the  vector  on  the  first  row  is  given  by  Marching  to 

the  last  row  for  each  unit  row  vector,  then  defines  the  influence  matrices  R21 
and  R22  for  the  last  two  rows  of  the  block.  If  error  vectors  C3,  64  and  £5 
represent  the  solution  error  on  the  second  and  last  two  rows  of  the  block, 
respectively,  then 

€4  =  €3  R21  (5.15) 

and 

65  =  €3  R22  (5.16) 

An  influence  matrix  $2  relating  the  solution  error  on  the  last  two  rows  of  the 
block  is  again  defined  so  that 

£5  =  £4  $2  (5.17) 

where  $2  =  ^21”^  ®22'  inverse  also  computed  using  82“^  =  ^22"^  ^21* 

For  the  last  block  (third  block  in  Fig.  2),  the  unit  vectors  are  marched 
to  the  next  to  last  row  and  residuals  computed,  defining  the  rows  of  a 
residual  matrix,  R31.  The  residual  matrix  then  relates  the  residual  solution 
error  p  on  the  last  interior  row  j  =  N-1  to  the  solution  error  £5  on  the 
second  row  j  =  IE(2)  of  the  last  block,  so  that 

P3  -  €5  ^31  (5.18) 
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5.3  The  Solution.  BSM2.  BSM3 

In  the  forward  sweep  an  approximate  solution  is  found  for  each  block, 
satisfying  the  boundary  conditions  and  an  initial  guess  field  on  the  last  row 
of  each  block.  This  particular  solution  satisfies  the  forcing  function  at  all 
points  except  on  the  last  Interior  row,  where  the  difference  equation 
satisfies  the  forcing  plus  a  residual. 

For  the  first  block,  we  want  to  obtain  the  particular  solution  which 
satisfies  an  assigned  initial  guess  value  on  the  interior  block  boundary 
j  =  IE(1).  The  boundary  values  and  the  initial  guess  on  the  second  row  are 
used  to  march  forward  to  the  block  boundary  and  the  error  CP3  from  the 
assigned  guess  value  on  the  block  boundary  is  computed.  The  required 
correction  on  the  second  row  is  then  given  by  eP^  =  CP3  Ri2"^.  applying  Eq. 
(5.13).  By  marching  the  error  forward,  the  corrected  particular  solution  is 
obtained.  The  process  is  repeated  for  several  iterations  to  reduce  the  round 
off  error. 

For  subsequent  interior  blocks,  the  particular  solution  on  the  first  two 
rows  of  the  block,  which  has  been  computed  in  the  previous  block,  is  similarly 
marched  forward  and  corrected  for  the  assigned  initial  guess  value  on  the  last 
row  of  this  block.  For  the  last  block  the  particular  solution  from  the 
previous  block  is  marched  to  the  last  interior  row  j  =  N-1  and  a  residual 
error  computed. 

In  the  backward  sweep,  error  vectors  are  computed  on  the  block  boundaries 
and  the  first  Interior  row  given  the  residual  error  on  the  last  interior  row, 
by  using  the  residual  and  influence  matrices  defined  for  the  blocks.  The 
solution  error  can  then  be  marched  forward  in  each  block  and  the  particular 
solution  corrected  to  give  the  solution. 

For  the  case  of  three  blocks,  as  shown  in  Fig.  3,  we  start  the  backward 
sweep  with  the  residual  error  p-^  given  on  the  last  interior  row  j  =  N-1.  The 
error  vectors  C4,  €5  on  the  first  two  rows  of  the  last  block  can  then  be 
computed  using  the  residual  matrix  R3X  and  the  influence  matrix  S2.  using 
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Fig.  3.  The  solution  errors  computed  in  the  backward  sweep  of  the  blocks, 
is  the  residual  error  on  the  row  j=N-l  and  c  represents  the  error  vectors  on 
the  rows  shown . 


65  =  ^3  R3i“^  (5.19) 

€4  =  £5  S2"^  (5.20) 

A  forward  marching  of  the  solution  error  then  corrects  the  solution  in  the 
last  block.  As  before,  several  iterations  are  used  to  reduce  round  off  error. 

In  the  backward  sweep  of  the  remaining  blocks,  the  particular  solution  is 
recomputed  on  the  last  row  of  each  block  and  the  solution  error  recomputed  on 
the  block  boundary.  For  the  second  block  in  Fig.  3  for  example,  the  error 
vector  €5  is  recomputed  by  backing  up  three  rows  into  the  second  block  and 
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repeating  a  small  segment  of  the  forward  sweep.  The  solution  error  is  then 
obtained  on  the  first  two  rows  of  the  block  using  Eqs.  (5.16)  and  (5.14), 

C3  »  €*5  8-22"^  (5.21) 

€2  “  €3  Si'l  (5.22) 

and  the  solution  corrected  throughout  the  block  as  above.  For  the  first  block 
the  recomputed  error  vector  e’3  on  the  last  row  provides  the  error  vector 
on  the  second  row  of  the  block  by  using  Eq  (5.13) 

€1  -  e’3  812“^  (5.23) 

and  the  solution  corrected  throughout  the  block. 

With  the  correction  of  the  solution  in  the  first  block  the  solution  is 
complete.  After  the  forward  and  backward  sweep,  the  remaining  residual  error 
accvunulates  on  the  interior  boundaries  of  the  blocks.  If  necessary,  the 
forward  and  backward  sweep  of  the  blocks  is  repeated  for  several  iterations  to 
reduce  any  remaining  error. 
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6.  DESCRIPTION  OF  TEST  EXAMPLE 


As  an  example  we  solve  Poisson’s  equation  ^  =  f  in  spherical  coordinates 

in  a  rectangular  domain  on  the  surface  of  a  spherical  earth.  We  define 
transformed  coordinates  x  and  y  by 
X  =  r  X  and  y  =  r  5 

where  r  »»  6371  km  is  the  average  radius  of  the  earth,  6  and  X  are  the  latitude 
and  longitude,  respectively  in  radians.  In  our  coordinates,  Poisson’s 
equation  is  written 


h  h 
X  y 


■  d 

h 

y 

8 

■ 

1 

.  0X 

.  h 

dx  . 

+ 

1 

.  h  8y  . 

j 

=  f(x,y),  a  b  (6.1) 

^a  ^  y  ^  ^b 


where  for  our  spherical  coordinates  the  map  factors  and  hy  are 
hjj  »  cos3  =  cos(y/r),  and  hy  =•  1. 

For  our  rectangular  domain  we  choose 

Xg  =  0  km,  x^  =  12,000  km,  y^  «  2,000  km  and  y^  =  7,000  km. 

We  want  to  force  a  wave  solution  on  a  zonal  background  field.  An  example  of 
such  a  solution  is  given  by 


^(x,y)  =  5.5  -  0.5  tanh(y  )  +  0.5  exp(-y  cos  ^  f  ^ 

P  P  ^ 


where 


(6.2) 


and  L  =  xfj  -  Xa  =  12,000  km,  y^.  =  4,500  km,  y^  =  1,500  km.  By  applying  Eq. 
(6.1).  the  appropriate  forcing  function  for  this  solution  is  then 


y  y 

f(x,y)  =  — ^  [  2  tanh(y  )  +  tan(  -  )  ]  sech^(y  )  + 

2  y  ^  ^ 

w 


2  ^w  y  2  2  w  X 

[  2  y  +  —  tan(  -  )  y  -  c  ]  exp(-y^)  cos(  — ) 

P  L  X  p  P  Lf 


(6.3) 


w 


where  c 


1  + 


2  2 

2  T  y 

•^w 


2  2  Y 

L  cos  (  ^  ) 
r 
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6.1  Dirichlet  Boundary  Conditions  in  v 


Our  boundary  conditions  are  taken  as  periodic  in  x  such  that 
^(x+Xj^,y)  =  ^(x+x^,y) 


and  in  this  first  case  we  treat  Dirichlet  boundary  conditions  at  the  y 
boundaries,  specifying  the  values  of  ^  at  y  *  y^^  and  y  “  y^.  That  is 


^(x,y  )  =  5.9655  +  0.0311  cos  ^  f- 

£t  M 

^(x.y^)  =  5.0344  +  0.0311  cos  -  " 


We  define  a  grid  of  M  =  82  by  N  =  51  points  (x^,  yj)  adding  an  extra  column 
of  points  at  x=xb+&x  for  the  cyclic  boundary  conditions.  The  grid  is  then 
Xi  =  (i-1)  Ax  +  x^  for  i  =  1,...,M 

yj  =  (j-1)  Ay  +  ya  for  j  =  1,...,N  (6.5) 

with  Ax  “  I  (H-Z)  =  150  km  and  Ay  ■  (yb-Xa^  /  ^ 

Our  Poisson’s  equation  (6.1)  can  now  be  written  in  finite  difference  form  as 


1 


h  .  h  . Ax 


xj  yj 


Vi  [  ^i+i.i  ~  ^i.i] 
h  .  Ax 


Ax 


+ 


h  .  h  . Ay 
xj  yj 


[  Vi+^  [  V.  i+1  '  V.  j]  Vi-4  [  V.i  '  V.1-l]  I 


h  .  ,  Ay 
yj+4 


h  .  ,  Ay 
yj-4  ^ 


(6.6) 


where  h^j+ij  =  (bxj+i  +  b^j )  /  2  and  with  the  boundary  conditions  (6.4) 
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2  T  X. 


^,1  - 

5.9655 

+  0.0311  cos 

1 

L 

2  T  X, 

^i,N 

5.0344 

+  0.0311  cos 

1 

L 

(6.7) 

and 

Vi.j 

^M.j 

The  difference  equation  (6.6)  can  further  be  written  as 


+ 


h  . 
xj 

Ax^ 

[  ^xi-t-%  [  ^i.i+l  "  ^i. 

i]  ^xi-^  [  ^i,.l  ■  ^i,i-i] 

h  . 

yj 

Ay^ 

i 

=  h  Ax^  f.  . 

XJ  l.J 

(6.7) 

Setting 

PHI(i.j) 

»  j  and  defining  the  coefficients  in  the  form 

for  the 

difference  equation 

(2.1)  we  have  for  i  =  2,...,M-1  and  j  =  2,...,N-1 

AX(i,j) 

= 

1. 

CX(i, j) 

= 

1. 

Ay(i,j) 

h^. 

(  h  .  +  h  .  ,) 

XI  xi-1' 

Ax^ 

h  . 

yj 

(  h  .  +  h  .  ,) 

yj  yj-1 

Ay^ 

CY(i, j) 

h^. 

(  h  +  h  . ) 

Xl+1  XT 

Ax^ 

(6.8) 

h  . 

yj 

(  h  .  ,  +  h  .) 
yj+1  y}' 

Ay^ 

BB(i,j) 

= 

-  2 

-  AY(i,j)  -  CY(i 

J) 

F(i.j) 

= 

h  ? 
XJ 

For  our  boundary  conditions  we  have  for  i  = 
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2  ir  X. 


PHI(i,l) 

s 

5.9655 

+  0.0311 

1 

cos  , 

with  All  = 

0.0 

i.. 

2  X 

(6.9) 

PHI(i,N) 

= 

5.0344 

+  0.0311 

cos  ^ 

with  AIN  ■= 

0.0 

and  ICY  = 

0, 

A21  = 

A2M  =  0 

.0  giving 

PHI(l,j) 

= 

PHI(M-1, 

.  j )  and 

PHI(M,j)  = 

PHI(2,j) 

(6.10) 

As  a  first 

guess  in  the 

interior 

of  the  domain,  we  initially  set  PHI 

to  the 

average  of 

the 

boundary  values  at 

y  =  Ya 

y  =  Yb*  giving 

PHI(i, j) 

= 

5.5  for 

i  =  2,  .. 

.  ,M-1  and  j 

=  2 . N-1. 

(6.11) 

With  a  call  to  SEVP,  the  solution  and  the  residual  error  are  computed.  The 
residual  error  gives  the  difference  between  the  forcing  and  the  discrete 
Laplacian  of  the  computed  solution.  The  maximum  of  the  absolute  values  of  the 
residual  errors,  divided  by  the  maximum  forcing  is  then  computed  and  printed. 
A  normalized  discretization  error  is  also  computed  by  taking  the  maximum  value 
of  the  absolute  differences  between  the  exact  analytic  solution  Eq.  (6.2)  and 
the  computed  solution  and  dividing  by  the  maximum  value  of  the  exact  solution. 
With  the  Cray  computer,  the  following  values  are  printed: 

Maximum  residual  error  (normalized  by  maximum  forcing)  *  5.9805  x  10"^^ 

Maximum  discretization  error  (normalized  by  maximum  value  of  exact  solution) 

-  8.3005  X  10-5. 


6 . 2  Neumann  Boundary  Conditions  in  y 

Our  boundary  conditions  are  again  taken  as  periodic  in  x  such  that 
^(x+Xj^,y)  =  ^(x+x^,y) 

but  in  this  case  we  treat  Neumann  boundary  conditions  at  the  y  boundaries, 
specifying  the  values  of  the  normal  gradient  of  ^  at  y  -  y^  and  y  *=  yb-  From 
our  definition  of  ^  the  y-gradient  of  ^  at  the  boundaries  y  *  ya  and  y  ■  yb  is 
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We  define  a  grid  of  M  =  82  by  N  =  53  points  (x^,  y j )  adding  an  extra  column 
of  points  at  x=X5+Ax  for  the  cyclic  boundary  conditions  and  two  rows  of  points 

at  y  =  ya  -  Ay/2  and  y  =  y^  +  Ay/2  for  the  Neumann  boundary  conditions.  The 

grid  is  then 

Xj^  =  (i-1)  Ax  +  Xa  for  i  = 

yj  =  (j-1)  Ay  +  ya  -  Ay/2  for  j  =  1,...,N  (6.13) 

with  Ax  =  (x^-Xa)  /  (M-2)  =  150  km  and  Ay  =  (yfYa)  /  (N-2)  =  100  km. 

Setting  PHI(i,j)  =  .  the  coefficients  for  the  difference  equation  are 

defined  by  Eq.  (6.8)  as  in  the  case  of  Dirichlet  boundary  conditions,  but  with 
the  map  factors  computed  for  the  new  grid.  For  our  boundary  conditions  we 
define 

Fll(i)  =  ■  ^ ■  sech^y 

2  y  ■'pa 

•^w  ^ 

FlN(i)  =  sech^y 

L  2  y  •'pb 

w 

with  All  =1.0  and  AIN  =  1.0 

and  ICY  =  0,  A21  =  A2M  =  0.0  giving 

PHI(l,j)  =  PHI(M-l.j)  and  PHI(M,j)  =  PHI(2,j)  (6.15) 

With  the  boundary  conditions  as  given,  the  solution  is  indeterminate.  To 
obtain  a  unique  solution,  we  require  a  value  for  the  solution  at  a  selected 
point  chosen  on  the  boundary  (preset  as  i  =  2,  j  =  1  in  the  SEVP  solver).  For 
the  most  efficient  convergence  of  the  solver  we  use  a  zero  value  for  PHI  at 
this  point  and  as  the  first  guess,  defining 


^pa 

exp( 

2  2  I-  x^ 

-V*  L  J 

Ay 

^w 

exp( 

2  2  »■ 
-ypb>  L  J 

Ay 
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PHI(i.j)  -  0.0  for  i  -  1 . M  and  j  =  1 . N.  (6.16) 

With  a  call  to  SEVP,  the  solution  and  the  residual  error  are  computed.  The 
residual  error  gives  the  difference  between  the  forcing  and  the  discrete 
Laplacian  of  the  computed  solution.  The  maximum  of  the  absolute  values  of  the 
residual  errors,  divided  by  the  maximum  forcing  is  then  computed  and  printed. 
A  normalized  discretization  error  is  also  computed  by  taking  the  maximum  value 
of  the  absolute  differences  between  the  exact  analytic  solution  and  the 
computed  solution  and  dividing  by  the  maximum  value  of  the  exact  solution. 
With  the  Cray  computer,  the  following  values  are  printed: 

Maximum  residual  error  (normalized  by  maximum  forcing)  =  8.0661  x  10“^^ 

Maximum  discretization  error  (normalized  by  maximum  value  of  exact  solution) 

*  8.8219  X  10-5. 
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APPENDIX:  LISTING  OF  TEST  PROGRAMS 


1.  Test  Example  with  Dirichlet  Boundary  Condition 

PROGRAM  DRIVRl 
C 

C  This  program  uses  SEVP  to  solve  Poisson’s  equation 
C  LAPLACIAN  OF  PHI  =  F 

C  with  cyclic  boundary  conditions  in  x  and  Dirichlet 
C  boundary  conditions  in  y. 

C 

PARAMETER  (M=82,  N=51) 

PARAMETER  (NBLK=(N+3) /9,  NBLK1=NBLK-1) 

PARAMETER  (M1=M-1,  M2=M-2,  N1=N-1,  N2=N-2) 

PARAMETER  (MN=M*N) 

DIMENSION  AX(M,N) . AY(M.N) ,CX(M,N) ,CY(M,N) ,BB(M,N) , 

1  RINV(M2,M2,NBLK) ,RINV1(M2,M2,NBLK1) ,RCOR(M,3) ,RTILDA(M2) , 

2  NBSIZ2(NBLK) ,IS(NBLK) ,IE(NBLK) ,SUMF(NBLK) ,IXNDX(M) , 

3  DUMMY(M2,M2) ,B1(M2) ,B2(M2) 

DIMENSION  F(M.N) .PHI(M.N) ,ERR(M,N) 

DIMENSION  Fll(M) ,F1N(M) , F21 (N) , F2M(N) 

DIMENSION  X(M) ,Y(N) ,YP1(N) ,YP2(N) 

DIMENSION  HXU(N) ,HXV(N1) ,HYU(N) ,HYV(N1) 

DIMENSION  AA(N) ,AC(N) ,CC(N) 

DIMENSION  T1(N) ,T2(N) ,T3(N) ,T4(N) 

DIMENSION  DUM3(M,N) 

C  SET  COEFFICIENTS  OF  DIFFERENCE  EQUATION  AND  BOUNDARY  CONDITIONS  TO 
C  ZERO  INITIALLY 

DATA  AX/MN*0.0/,CX/MN*0.0/,AY/MN*0.0/,CX/MN*0.0/,BB/MN*0.0/, 

1  F/MN*0.0/,F11/M*0.0/,F1N/M*0.0/,F21/N*0.0/,F2M/N*0.0/ 

C 

PI  =  A.0*ATAN(1.0) 

AR  =  6371. 

XL  =  12000.0 
YL  =  5000.0 

DEFINE  GRID 

DELX  =  XL / FLOAT (M- 2) 

DELY  =  YL / FLOAT (N-1) 

PRINT  9 90, DELX, DELY 

990  FORMATdH  ,  5HDELX=,  1PE13 . 5, 5X,  5HDELY=,  1PE13 . 5) 

DO  5  J=1,N 

Y(J)=DELY*FLOAT(J-1)+2000. 

5  CONTINUE 
DO  6  1=1, M 
X ( I ) =DELX*FLOAT ( I- 1 ) 

6  CONTINUE 
DELXSQ=DELX*DELX 
DELYSQ=DELY*DELY 

C  DEFINE  MAP  FACTORS 
DO  10  J=1,N 
HXU(J)=COS(Y(J)/AR) 

HYU(J)=1.0 
10  CONTINUE 
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DO  20  J-1,N1 

HXV(J)-0.5*(HXU(J)+HXU(J+1)) 

HYV(J)-1.0 
20  CONTINUE 

DEFINE  FORCING  FUNCTION 
XC=  0 . 0 
YC  =  4500.0 
YW  =  1500.0 
C1=YW/AR 
C2=1./(YW*YW) 

C3=2.*PI*PI*(YW/XL)*(YW/XL) 

DO  60  J=1,N 
YP1(J)=(Y(J)-YC)/YW 
YP2(J)  =  Y(J)/AR 

60  CONTINUE 

DO  61  J=1.N 

T1(J)  =  TANH(YP1(J)) 

T2(J)  -  EXP(-YP1(J)*YP1(J)) 

61  CONTINUE 

DO  62  J=1,N 

CSQ=COS ( YP2 ( J ) ) *C0S ( YP2 ( J ) ) 
T3(J)=2.*T1(J)+C1*TAN(YP2(J)) 

T4(J)=  2.*YP1(J)*YP1(J)  -  1.0  +  C1*YP1(J)*TAN(YP2(J)) 

1  -  C3/CSQ 

62  CONTINUE 

FORCING  FUNCTION  DEFINED  AT  INTERIOR  POINTS  OF  DOMAIN  ONLY 
DO  65  J=2,N-1 
DO  65  1=2, M-1 

F(I,J)  =  0.5*C2*T3(J)*(1.-T1(J)*T1(J)) 

1  +  C2*T4(J)*T2(J)*COS(2.*PI*(X(I)-XC) /XL) 

65  CONTINUE 


DIRICHLET  BOUNDARY  CONDITIONS  IN  Y 
All  =0.0 
AIN  =0.0 
DO  100  1=1, M 

PHI(I,1)=5.5-0.5*T1(1)+0.5*T2(1)*COS(2.*PI*(X(I)-XC)/XL) 
PHI(I,N)=5.5-0.5*T1(N)+0.5*T2(N)*COS(2.*PI*(X(I)-XC)/XL) 
100  CONTINUE 

FOR  DIRICHLET  B.C.  INITIALLY  SET  PHI  IN  INTERIOR  TO  AVERAGE 
OF  BOUNDARY  VALUES 
AVPH  =0.0 
DO  110  1=2, Ml 

AVPH  =  AVPH+PHI(I,1)+PHI(I,N) 

110  CONTINUE 

A VPH=AVPH / FLOAT ( 2  * ( M- 2 ) ) 

DO  120  J=2,N1 
DO  120  1=2, Ml 
PHI(I,J)»  AVPH 
120  CONTINUE 
C  FOR  CYCLIC  BOUNDARY  CONDITIONS 
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ICYC  -  1 
A21  -  0.0 
A2M  -0.0 
DO  130  J=1,N 
PHKl.J)  =  PHI(M-1,J) 

PHI(M,J)  -  PHI(2,J) 

130  CONTINUE 
C 

C  DEFINE  COEFFICIENTS  ON  L.H.S.  OF  DIFFERENCE  EQUATION. 

C  COEFFICIENTS  DEFINED  IN  INTERIOR  OF  DOMAIN  ONLY. 

DO  25  J=2,N-1 

AA(J)=( (HXU( J) *HXV( J-1) ) / (HYV(J-1)*HYU( J) ) )*(DELXSQ/DELYSQ) 

CC(J)=((HXU(J)*HXV(J))/(HYV(J)*HYU(J)))*(DELXSQ/DELYSQ) 

AC(J)^-AA(J)-CC(J) 

25  CONTINUE 

DO  165  I  =  2,M-1 
DO  165  J  =  2,N-1 
AX(I.J)  =  1.0 
CX(I,J)  =  1.0 
AY(I,J)  =  AA(J) 

CY(I,J)  =  CC(J) 

BB(I,J)  =  AC(J)-2.0 
165  CONTINUE 

C  DEFINE  R.H.S.  OF  DIFFERENCE  EQUATION,  SCALING  FORCING  FUNCTION 
DO  170  J=2,N-1 
DO  170  1=2, M-1 
CON  =  HXU(J)*HXU(J)*DELXSQ 
F(I,J)  =  F(I,J)*CON 
170  CONTINUE 
FMAX  =  0. 

DO  305  J=2,N-1 
DO  305  1=2, M-1 

FMAX  =  AMAX1(FMAX,ABS(F(I,J))) 

305  CONTINUE 


TOLERANCE  FOR  SOLVER 
TOL  -  l.OE-8 

IFLG  =  1 

CALL  SEVP ( AX , AY , CX , CY , BB , RINV , RINVl , RCOR, NBSIZ2 , IS , IE , IXNDX , 
1SUMF,DUMMY,RTILDA,PHI,F,ERR,F11,F1N,F21,F2M,M,N,M1,N1,M2,N2, 
2NBLK , NBLKl , B1 , B2 , All , AIN , A21 , A2M, ICYC , TOL , IFLG ) 

PRINT  710 
PRINT  711,  NBSIZ2 

710  FORMAT (IX, ’BLOCK  SIZE  GIVEN  BY  NBSIZ2(l)+2  (1ST  BLOCK)  OR’, 

1  ’  NBSIZ2()+3  (OTHER  BLOCKS)’/, IX, ’  NBSIZ2:  ’) 

711  FORMAT  (IX,  1014) 

RESIDUAL  ERROR  OF  LAPLACIAN 
ERMAX=0 . 0 
DO  750  1=2, Ml 
DO  750  J=2,N1 
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ERMAX  -  AMAX1(ERMAX,ABS(ERR(I,J)/FMAX)) 

750  CONTINUE 

PRINT  751,  ERMAX 

751  PORMATdX,’  MAXIMUM  RESIDUAL  ERROR  [(FORCING  -  LAPLACIAN) /MAX’ , 

1  ’  FORCING  ]  »  ’ , 1PE12 . A ) 

C 

C  DISCRETIZATION  ERROR 
PMAX  =  0. 

DO  3000  J=1,N 
DO  3000  1=1, M 
PHIX=  5.5-0.5*Tl(J) 

1  +  0.5*T2(J)*C0S(2.*P1*(X(I)-XC)/XL) 

PMAX  =  AMAX1(PMAX,ABS(PHIX)) 

DUM3(I, J)=PHI(I, J)-PHIX 
3000  CONTINUE 
DSMAX=0 . 0 
DO  3007  1=2, Ml 
DO  3007  J=2,N1 

DSMAX  =  AMAX1(DSMAX,ABS(DUM3(I,J)/PMAX) ) 

3007  CONTINUE 

PRINT  732,  DSMAX 

732  FORMAT ( IX, ’  MAXIMUM  DISCRETIZATION  ERROR  [(COMPUTED  -  ANALYTIC’, 
1  ’  PHI) /MAX  ANALYTIC  PHI]  =  ’,1PE12.A) 

C 

600  CONTINUE 
END 
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2.  Test  Example  with  Neumann  Boundary  Condition 

PROGRAM  DRIVR2 
C 

C  This  program  uses  SEVP  to  solve  Poisson’s  equation 
C  LAPLACIAN  OF  PHI  =  F 

C  with  Neumann  boundary  conditions  in  y 
C  and  Cyclic  boundary  conditions  in  x 
C 

PARAMETER  (M=82,  N=53) 

PARAMETER  (NBLK-l+(N-5+3) /9.  NBLK1=NBLK- 1 ) 

PARAMETER  (M1=M-1,  M2=M-2,  N1=N-1.  N2-N-2) 

PARAMETER  (MN=M*N) 

DIMENSION  AX(M,N) ,AY(M,N) ,CX(M,N) .CY(M,N) ,BB(M,N) , 

1  RINV (M2 , M2 , NBLK) . RINVl (M2 ,M2 , NBLKl ) , RCOR (M, 3 ) , RTILDA (M2 ) . 

2  NBSIZ2 (NBLK) .IS (NBLK) , IE (NBLK) ,SUMF (NBLK) ,IXNDX(M) , 

3  DUMMY(M2.M2)  ,B1(M2)  ,B2(M2) 

DIMENSION  F(M,N) ,PHI(M,N) ,ERR(M,N) 

DIMENSION  Fll (M) , FIN (M) , F21 (N) , F2M(N) 

DIMENSION  X(M) ,Y(N) ,YP1(N) ,YP2(N) 

DIMENSION  HXU(N) ,HXV(N1) ,HYU(N) ,HYV(N1) 

DIMENSION  AA(N) .AC(N) ,CC(N) 

DIMENSION  T1(N).T2(N),T3(N),T4(N) 

DIMENSION  DUM3(M,N) 

C  SET  COEFFICIENTS  OF  DIFFERENCE  EQUATION  AND  BOUNDARY  CONDITIONS 
C  TO  ZERO  INITIALLY. 

DATA  AX/MN*0 . 0 / , CX/MN*0 . 0 / , AY/MN*0 . 0 / , CX/MN*0 . 0 / , BB/MN*0 . 0 / , 
1  F/MN*0 .01, Fll /M*0 .01, F1N/M*0 . 0 / , F21 /N*0 . 0 / , F2M/N*0 . 0 / 

C 

PI  =  4.0*ATAN(1.0) 

AR  =  6371. 

XL  =  12000.0 
YL  =  5000.0 

DEFINE  GRID 

DELX  =  XL / FLOAT (M- 2) 

DELY  =  YL/FLOAT(N-3) 

PRINT  990, DELX, DELY 

990  FORMATdH  ,  5HDELX=,  1PE13 . 5 , 5X,  5HDELY-,  1PE13 . 5) 

DO  5  J=1,N 

Y(J)=DELY*(FL0AT(J-l)-0.5)  +  2000.0 

5  CONTINUE 
DO  6  1=1, M 
X ( I ) =DELX*FLOAT ( I - 1 ) 

6  CONTINUE 
DELXSQ=DELX*DELX 
DELYSQ=DELY*DELY 

DEFINE  MAP  FACTORS 
DO  10  J=1,N 
HXU(J)=C0S(Y(J)/AR) 

HYU(J)=1.0 
10  CONTINUE 
DO  20  J=1,N1 
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HXV(J)*0.5*(HXU(J)+HXU(J+1)) 

HYV(J)=1.0 
20  CONTINUE 

DEFINE  FORCING  FUNCTION 
XC=  0.0 
YC  =  4500.0 
YW  =  1500.0 
Cl=YVr/AR 
C2=l. / (YW*YW) 

C3=2.*PI*PI*(YW/XL)*(YW/XL) 

DO  60  J=1,N 
YP1(J)=(Y(J)-YC) /YW 
YP2(J)  =  Y(J)/AR 

60  CONTINUE 

DO  61  J=1,N 

T1(J)  =  TANH(YP1(J)) 

T2(J)  =  EXP(-YP1(J)*YP1(J) ) 

61  CONTINUE 

DO  62  J=1,N 

CSQ=COS(YP2(J) )*C0S(YP2(J) ) 
T3(J)=2.*T1(J)+C1*TAN(YP2(J)) 

T4(J)=  2.*YP1(J)*YP1(J)  -  1.0  +  C1*YP1(J)*TAN(YP2(J)) 
1  -  C3/CSQ 

62  CONTINUE 

FORCING  FUNCTION  ONLY  DEFINED  IN  INTERIOR  OF  DOMAIN. 

DO  65  J=2,N-1 
DO  65  1=2, M-1 

F(I,J)  =  0.5*C2*T3(J)*(1.-T1(J)*T1(J)) 

1  +  C2*T4(J)*T2(J)*COS(2.*PI*(X(I)-XC) /XL) 

65  CONTINUE 


NEUMANN  BOUNDARY  CONDITIONS  IN  Y 
All  =1.0 
AIN  =1.0 

Yll  =  (Y(l)+(DELY/2.0)-YC) /YW 
YIN  =  (Y(N)-(DELY/2.0)-YC)/YW 
Til  =  TANH(Yll) 

TIN  =  TANK (YIN) 

T21  =  EXP(-Y11*Y11) 

T2N  =  EXP(-Y1N*Y1N) 

DO  150  I  =  2, M-1 

Fll(I)  =  -  (DELY/YW)*(  0.5*(1.0-T11*T11) 

1  +Y11*T21*C0S(2.*PI*(X(I)-XC) /XL)  ) 

FIN(I)  =  -  (DELY/YW)*(  0 . 5* (1 . 0-TlN*TlN) 

1  +Y1N*T2N*C0S(2.*PI*(X(I)-XC)/XL)  ) 

150  CONTINUE 
C  CYCLIC  BOUNDARY  CONDITIONS  IN  X 
ICYC  =  1 
A21  =  0.0 
A2M  =0.0 
C 

C  FOR  INDETERMINATE  NEUMANN  BOUNDARY  CONDITIONS,  A  BOUNDARY  VALUE 
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C  MUST  BE  PROVIDED  AT  A  SINGLE  BOUNDARY  POINT  (IX.JX)  FOR  THE  SOLVER  TO 
C  OBTAIN  A  SOLUTION.  THE  BOUNDARY  POINT  SELECTED  IN  SEVP  IS  IX*2,  JX-1. 
C  FOR  THE  EFFICIENT  CONVERGENCE  OF  THE  SOLUTION,  WE  USE  A  ZERO  VALUE 
C  FOR  THE  SINGLE  BOUNDARY  POINT  AND  AS  THE  FIRST  GUESS. 

C 

DO  120  J=1,N 
DO  120  I=-1,M 
PHI(I,J)-  0.0 
120  CONTINUE 
C 

C  DEFINE  COEFFICIENTS  ON  L.H.S.  OF  DIFFERENCE  EQUATION. 

C  COEFFICIENTS  DEFINED  IN  INTERIOR  OF  DOMAIN  ONLY. 

DO  25  J»2,N-1 

AA(J)*( (HXU(J)*HXV(J-1) ) / (HYV(J-1)*HYU(J) ) ) * (DELXSQ/DELYSQ) 

CC(J)-t (HXU(J)*HXV(J))/ (HYV(J)*HYU{J)))*(DELXSQ/DELYSQ) 
AC(J)=-AA(J)-CC(J) 

25  CONTINUE 

DO  165  I  -  2.M-1 
DO  165  J  -  2.N-1 
AX(I,J)  -  1.0 
CX(I.T)  =  1.0 
AY(I.J)  =  AA(J) 

CY(I.J)  »  CC(J) 

BB(I,J)  -  AC(J)-2.0 
165  CONTINUE 

C  DEFINE  R.H.S.  OF  DIFFERENCE  EQUATION,  SCALING  FORCING 
DO  170  J=2,N-1 
DO  170  1-2, M-1 
CON  »  HXU(J)*HXU(J)*DELXSQ 
F(1,J)  -  F(I,J)*CON 
170  CONTINUE 
FMAX  -  0. 

DO  305  J-2,N-1 
DO  305  1=2, M-1 

FMAX  =  AMAX1(FMAX,ABS(F(I,J))) 

305  CONTINUE 


TOLERANCE  FOR  SOLVER 
TOL  =  l.OE-8 

IFLG  =  1 

CALL  SEVP ( AX , AY , CX , CY , BB , RINV , RINVl , RCOR ,NBSIZ2,IS,IE, IXNDX , 
1SUMF,DUMMY,RTILDA,PHI,F,ERR,F11,F1N,F21,F2M,M,N,M1,N1,M2,N2, 
2NBLK , NBLKl , B1 , B2 , All , AIN , A2 1 , A2M, ICYC , TOL , IFLG ) 

PRINT  710 

710  FORMAT (IX, ’BLOCK  SIZE  GIVEN  BY  NBSIZ2(l)+2  (1ST  BLOCK)  OR’, 

1  ’  NBSIZ2()+3  (OTHER  BLOCKS) ’/, IX, ’  NBSIZ2;  ’) 

PRINT  711,  NBSIZ2 

711  FORMAT (IX, 1014) 

RESIDUAL  ERROR  OF  LAPLACIAN 
FNORM-FMAX 
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ERMAX*0.0 
DO  750  1=2, Ml 
DO  750  J=2,N1 

ERMAX  =  AMAXl{ERMAX,ABS(ERR(I,J)/FNORM)) 

750  CONTINUE 

PRINT  751.  ERMAX 

751  FORMATCIX,’  MAXIMUM  RESIDUAL  ERROR  [(FORCING  -  LAPLAC IAN) /MAX’ , 

1  ’  FORCING  ]  =  ■,1PE12.4) 

C 

C  DIFFERENCE  AT  PRESCRIBED  BOUNDARY  POINT  (IX.JX) 

PHI21=5.5-0.5*T1(1)+0.5*T2(1)*COS(2.*PI*(X(2)-XC)/XL) 

PDIFF  =  PHI21  -  PHI(2,1) 

C  DISCRETIZATION  ERROR 
PMAX  =  0. 

DO  3000  J=1,N 
DO  3000  1=1, M 
PHIX=  5.5-0.5*Tl(J) 

1  +  0.5*T2(J)*COS(2.*PI*(X(I)-XC)/XL) 

DUM3 ( I , J ) =PHI ( I . J ) -PHIX+PDIFF 
PMAX  =  AMAX1(PMAX.ABS(PHIX)) 

3000  CONTINUE 

PNORM  =  PMAX 
DSMAX=0.0 
DO  3007  1=2, Ml 
DO  3007  J=2,N1 

DSMAX  =  AMAX1(DSMAX,ABS(DUM3(I,J)/PN0RM)) 

3007  CONTINUE 

PRINT  732,  DSMAX 

732  FORMATCIX,’  MAXIMUM  DISCRETIZATION  ERROR  [(COMPUTED  -  ANALYTIC’, 
1  ’  PHI) /MAX  ANALYTIC  PHI]  =  ',1PE12.4) 

C 

600  CONTINUE 
END 
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